/^DA  01  802  0 


implementing  VEHICl-E  ROUTING  ALGORITHMS 

by 

BRUCE  L.  GOLDEN 
THOMAS  L.  MAGNANTI 
and 

HIEN  Q.  NGUYEN 


leGhniedl  Report  No,  115 
OPERATIONS  RESEARCH  CENTER 


-^10  D C'V 


I^A^SACMIljSEm  iNSTITUTl 
^ OF 

tECHNOfcQGYr’i^^M.S^S 


September  1Qf5 


Ji-  \ ^ /( 


DISCLAIMER  NOTICE 


THIS  DOCUMENT  IS  BEST  QUALITY 
PRACTICABLE.  THE  COPY  FURNISHED 
TO  DTIC  CONTAINED  A SIGNIFICANT 
NUMBER  OF  PAGES  WHICH  DO  NOT 
REPRODUCE  LEGIBLY. 


■Ball 


BY  - TTT- 

. BlSYRlB«Tl«S/Ay/.IUclUn'  cora 


□<a 


Unclassified 

SECURITY  CLASSIFICATION  OF  THIS  PACE  fW»i«n  Dels  Entered) 

REPORT  DOCUMENTATION  PAGE  before'^completS 

I.  REPORT  NUMBER  |2.  GOVT  ACCESSION  NO.  3.  RECIPIENT’S  CAT  ALOG  NUMBER 


I.  REPORT  NUMBER 


Technical  Report  No.  115 

*.  title  fBnd  Si  — T-VPF  fiF  RFPnPT  A PFtiirin  rn\jtrarn 

/ / Technical  ^ep^t  I 

L j IMPLEMENTING  VEHICLE  ROUTING  ALGORITHMS^  ) L/ ^^—fhrpt-crohTrr^r^TT^^ 

f ^ ^ ^ J 

' ) R PERFORMING  ORG.  REPORT  NUMBER 


l_.  AUTHORfj) 


ID  '(  Bruce  L.  /Goldc’  >,  \ /,/]  - , ^ 

y'  Thomas  L./Magnanti  y •^7  (7 

Hien  Q./Nguven  y [ 

3 PERFORMING  ORG ANIzllTION  NAME  ANO  ADDRESS 

M.I.T.  Operations  Research  Center  ^ 

77  Massachu.setts  Avenue,  Room  24-215  [1 1 

Cambridge,  M\  02139  ^ 

II.  CONTROLLING  O FICE  NAME  ANO  ADDRESS 

O.R.  Branch.  ONR  Navy  Dept. 

800  North  Quincy  Street  I 

Arlington,  VA  22217  4 

t4  MONITORING  AGENCY  NAME  6 AOORESSf// di//«rcfif  from  Omffoffin/}  0//icoy 


B.  CONTRACT  OR  GRANT  NUMBERf*; 

3 NQ0D14-75-C-^556./  

fDAHCdQ~03'C- 

CS^iMTO5RiSr5555^^ 

AREA  B WORK  UNIT  NUMBERS 

^ NR-347-ji(27  j 

12.  REPORT  DATE 

) SepfeSSi^STr 

13.  NUMBER  OF  pages 

46 

IS.  SECURITY  CLASS,  fo/ «i/»  reporO 

Unclassified 


ISn.  DECLASSIFICATION/ DOWN  GRADING 
SCHEDULE 


1 16.  Distribution  statement  (ot  ihie  Report) 


^-pprc^od  for  pubuHr 


17  Distribution  statement  (at  the  abetmcl  entered  in  Block  20,  It  dlllerenl  Itom  Report) 


18.  supplementary  notes 


19.  KEY  WORDS  fConl/nuo  on  reverse  side  tt  nocosooiy  and  Identlly  by  block  number) 


Heuristic  Programming 
Vehicle  Routing 
Multi-Depot  Routing  Algorithm 


ZO.j^BSTRACT  (Continue  on  reverae  aide  If  neceaamry  end  Identity  by  block  number) 

Heuristic  programming  algorithms  frequently  address  large  problems  and  require 
manipulation  and  operation  on  massive  data  sets.  The  algorithms  can  be  im- 
proved by  using  efficient  data  structures.  With  this  in  mind,  we  consider 
heuristic  algorithms  for  vehicle  routing,  comparing  techniques  of  Clarke  and 
Wright,  Gillett  and  Miller,  and  Tyagi,  and  presenting  modifications  and  ex- 
tensions which  permit  problems  involving  hundreds  of  demand  points  to  be 
solved  in  a matter  of  seconds.  In  addition,  a multi-depot  routing  algorithm — ' 


DD  4 


FORM 
JAN  73 


EDITION  OF  I NOV  65  IS  OBSOLETE 


Unclassified  0 

SECURITY  CLASSIFICATION  ©"f  THIS  PAGE  flWion  Date  Entered)/ 


■ '■ 

jUnclassifieJ 


SECURITY^LASSIFICATION  of  this  PAQCfWiwi  Dmtm 


(9^  is  developed.  The  results  are  illustrated  with  a routing  study  for  an 
urban  newspaper  with  an  evening  circulation  exceeding  100,000. 


Unclassified 


SECURITY  CLASSIFICATION  OF  TH'S  PAGEriWion  Doik  Bnleted) 


IMPLEMENTING  VEHICLE  ROUTING  ALGORITHMS 


by 

BRUCE  L.  GOLDEN* 
THOMAS  L.  MAGNANTI** 
and 

HIEN  Q.  NGUYEN 


Technical  Report  No.  115 


Operations  Research  Center 
Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts  02139 


September  1975 


niibUc  relcoco; 


Supported  in  part  by  the  American  Newspapers  Publishers'  Association 
and  the  Department  of  Transportation  under  Contract  DOT-TSC-1058. 

Supported  in  part  by  the  Office  of  Naval  Research  under  Contract 
NQ0014-75-C-0056  and  in  part  by  the  Army  Research  Office  under 
CorvEHcF  TDimeOA- 7 3-C-003  2 . 


ABSTRACT 


Heuristic  programning  algorithms  frequently  address  large  prob- 
lems and  require  manipulation  and  operation  on  massive  data  sets.  The 
algorithms  can  be  improved  by  using  efficient  data  structures.  With  this 
in  mind,  we  consider  heuristic  algorithms  for  vehicle  routing,  comparing 
techniques  of  Clarke  and  Wright,  Gillett  and  Miller,  and  Tyagi,  and  pre- 
senting modifications  and  extensions  which  perndt  problems  involving 
hundreds  of  demand  points  to  be  solved  in  a matter  of  seconds.  In  ad- 
dition, a multi-depot  routing  algorithm  is  developed.  The  results  are 
illustrated  with  a routing  study  for  an  urban  newspaper  with  an  evening 
circulation  exceeding  100,000. 
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I . INTRODUCTION 


An  essential  element  of  any  logistics  system  is  the  allocation  and 
routing  of  vehicles  for  the  purpose  of  collecting  and  delivering  goods  and 
services  on  a regular  basis.  Common  examples  include  newspaper  delivery  [29], 
schoolbus  routing  [8],  municipal  waste,  collection  [7],  fuel  oil  delivery  [25], 
and  truck  dispatching  in  any  of  a number  of  industries.  The  system  may  in- 
volve a single  depot  or  multiple  depots;  the  objectives  may  be  aimed  at 
cost  minimization  (distribution  costs,  and  vehicle  or  depot  acquisition  costs) 
or  service  improvement  (increasing  distribution  capacities,  ’•educing  distri- 
bution time,  and  related  network  design  issues).  Constraints  may  be  imposed 
upon 

(i) 

(ii) 

(iii) 

(iv) 

(v) 

and  (vi) 

Applications  emphasizing  various  of  these  characteristics  lead  to  a 
continuum  of  overlapping  models  including  location  models  such  as  median 
and  minimax  location,  scheduling  models  such  as  crew  scheduling,  distribution 
models  such  as  minimum  cost  flow  or  shortest  paths,  or  location-distribution 


the  depots  (numbers,  possible  locations,  and  production 
capabilities) , 

the  vehicle  fleet  (types  and  numbers  of  vehicles,  and 
vehicle  capacities), 

the  delivery  points  (demand  requirements,  service  con- 
straints on  delivery  time,  and  order  splitting), 

the  routing  structure  (maximum  route  time  or  route 
distance,  link  capacities,  and  preferences  for  radial 
routes,  peripheral  routes,  or  routes  with  points  closer 
together) , 

operator  scheduling  and  assignments  (union  regulations) , 

system  dynamics  (inventory  holdings,  and  distribution 
or  acquisition  lag  times) . 
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models  such  as  warehouse  location.  In  this  paper,  we  consider  models  of  the 
routing  type  with  fixed  depots,  constrained  for  the  most  part  by  fleet, 
delivery  point,  and  route  structure  restrictions.  Our  purposes  are  three-fold. 
By  synthesizing  and  extending  the  results  presented  in  Golden  [29]  and 
Nguyen  [54],  we  first  formulate  integer  programming  models  of  these  problems, 
building  upon  the  traveling  salesman  problem  as  the  "core"  model.  We  then 
review  the  most  promising  heuristic  solution  procedures  for  these  problems 
suggested  in  the  literature.  Finally,  we  suggest  particularly  fast  implemen- 
tations of  the  Clarke-Wrlght  heuristic  procedure.  Our  method  uses  efficient 
data  handling  capabilities  such  as  heap  structures  to  enhance  both  computations 
and  storage.  Our  results  indicate  that  algorithmic  modifications  of  this 
nature  can  greatly  extend  the  size  of  problems  amenable  to  analysis. 

We  have  applied  the  algorithm  to  a distribution  system  for  an  urban 
newspaper  with  an  evening  circulation  exceeding  100,000.  This  problem 
contains  nearly  600  drop  points  for  newspaper  bundles  and  was  solved  with  20 
seconds  of  execution  time  on  an  IBM  370/168.  Fast  computations  of  this 
nature  permit  the  algorithm  to  be  used  as  an  operational  tool  and  to  be  used 
interactively  to  study  modeling  assumptions  such  as  demand  patterns.  It 
could  be  used,  for  example,  as  in  our  newspaper  study,  to  investigate  perfor- 
mance measures  associated  with  operating  policies  like  the  pre-routing  of 
important  drop  points  on  special  routes.  In  addition,  past  implementation 
suggests  that  the  heuristic  might  be  used  as  a subroutine  for  more  general 
problems  in  order  to  incorporate  other  of  the  distribution  characteristics 
delineated  above,  or  to  perform  on-line  dynamic  routing. 

We  believe  that  our  computational  experience  is  indicative  of  implementa- 
tions for  heuristic  algorithms  in  general.  Because  these  algorithms  frequently 
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replace  computationally  expensive  or  intractable  formal  optimization 
procedures  by  simple,  but  usually  extensive,  addition  and  comparison  operations, 
data  organization  can  be  important.  Improving  data  access  for  each  basic 
operation  can  lead  to  order  of  magnitude  improvements  in  overall  running  time. 


II,  FORMULATIONS 


The  Traveling  Salesman  Problem 

Most  vehicle  routing  models  are  variants  and  extensions  of  the  ubiquitous 
Traveling  Salesman  Problem;  for  a very  thorough  overview  of  this  problem  see 
Bellmore  and  Nemhauser  [6].  Suppose  we  are  given  the  matrix  of  pairwise 
distances  or  costs  d^^  between  node  i and  node  j for  the  n nodes  1,  2,  ...,  n. 
We  assume  d^^  = " for  i = 1,  2,  n.  The  problem  is  to  form  a tour  of  the 

n nodes  beginning  and  ending  at  the  origin,  node  1,  which  gives  the  minimum 
total  distance  or  cost.  For  notation,  let 

n = the  number  of  nodes  in  the  network 
V = the  set  of  nodes 


if  arc  (i,  j)  is  in  the  tour 
otherwise 


d^j  = the  distance  on  arc  (i,  j). 

An  assignment-based  ‘'ormulation  of  the  problem  selects  the  matrix  X = (x_  ) of 
decision  variables  solving: 


n n 


Minimize 


I Z d..  X,  . 
i=i  j-i  ij 


(1.1) 
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subjccC  Co 


(j  = 1,  ....  n) 


(1.2) 


Z = 1 

j=l  ^ 


.X  = (x^j)  £ S 


(1.4) 


X . . = 0 or  1 


(i  Ij  ....  nj 

j = 1,  ....  n) 


(1.5) 


The  set  S is  selected  to  prohibit  subtour  solutions  satisfying  the  assign- 
ment constraints  (1.2),  (1.3),  and  (1.5).  Several  alternates  have  been  pro- 
posed for  S including 

(1)  S = {(Xj^.):  Z Z x^.  1 for  every  nonempty 

^ ieQ  j«!Q  proper  subset  Q of  V } ; 


(2)  S = 


{(x^J:  Z Z 

ieq  j£Q 


£ |q|  - 1 for  every  nonempty 

subset  Q of 


5CL 

. 3,  ....  ’ 


(3)  S = {(x^J:  y - y_  + nx  < n-1  for  2 _<  i ^ j < n 


1 J 


for  some  real  numbers  y./, 

1 


Note  that  S contains  nearly  2 subtour  breaking  constraints  in  (1)  and  (2), 

2 

but  only  n - 3n  + 2 constraints  in  the  ingenious  formulation  (3)  proposed  by 
Miller,  Tucker,  and  Zemlin  [51].  Algorithmically,  the  constraints  in  (2) 
have  proved  very  useful,  however,  for  Lagrangian  approaches  to  the  Traveling 
Salesman  Problem,  as  initiated  by  Held  and  Karp  [32],  [33].  For  other  approaches, 
both  heuristic  and  optimal,  see  Lin  [46],  Lin  and  Kernighan  [47],  Christofides 
and  Eilon  [16],  Little  [48],  and  Shapiro  [64]. 


I 


-5- 


The  Traveling  Salesman  Problem  can  be  interpreted  as  a vehicle  routing 
model  with  one  depot  and  with  one  vehicle  whose  capacity  exceeds  total  demand. 
This  model  can  be  extended  by  considering  more  vehicles,  more  depots,  different 
vehicle  capacities,  and  additional  route  restrictions. 

The  Multiple  Traveling  Salesmen  Problem 

The  Multiple  Traveling  Salesmen  Problem  (MTSP)  is  a generalization  of 
the  Traveling  Salesman  Problem  (TSP)  and  comes  closer  to  accomodating  more 
real-world  problems;  here  there  is  a need  to  account  for  more  than  one  salesman 
(vehicle) , Multiple  Traveling  Salesmen  Problems  arise  in  many  sorts  of 
scheduling  and  sequencing  applications.  For  example,  the  framework  could  be 
used  to  develop  the  basic  route  structure  for  a pickup  or  delivery  service 
(perhaps  a schoolbus  or  rural  bus  service) ; it  has  proved  to  be  an  appropriate 
model  for  the  problem  of  bank  messeng  r scheduling,  where  a crew  of  messengers 
picks  up  deposits  at  branch  banks  and  returns  them  to  the  central  office  for 
processing  [65]. 

Given  M salesmen  and  n nodes  in  a network  the  MTSP  is  to  find  M subtours 
(each  of  which  includes  the  origin)  such  that  every  node  (except  origin)  is 
visited  exactly  once  by  exactly  one  s.Tlesman,  so  that  the  total  distance  travel- 
ed by  all  M salesmen  is  minimum.  A MTSP  formulation  is  displayed  below. 


: Minimize 

e" 

d^^ 

X,  . 

1 

i=l 

j=l 

ij 

ij 

r.n 

1 

fM 

if 

3=^1 

j “2 y 3 y •••}  n 

subject  to 

•i 

i 

E 

i=l 

^j 

i 

if 

e" 

V 

1 

^M 

if 

3=1 

I 

Lt 

j=l 

A , , 

13 

" ^i 

1 

1 

1 

if 

1^2  y 3 y * « * y I\ 

I 


!) 


X 
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X =■  (x  ) £ S (2.4) 

Ij 

Xj,.  = 0 or  1 (i  = 1,  n;  (2.5) 

j ~ ...»  n) 

for  any  choice  of  the  set  S that  breaks  subtours  which  do  not  include  the 
origin.  In  particular,  any  of  the  three  choices  given  previously  for  S can  be 
used . 

Svestka  and  Huckfeldt  [65]  present  a Hiller-Tucker-Zemlin-like  formulation 
for  S and  apply  a subtour  elimination  type  branch  and  bound  procedure  using 
Bellmore  and  Malone  branching  to  obtain  the  optimal  solution;  mean  run  time  for 
55  city  problems  is  one  minute.  Three  different  papers  published  in  1973 
and  1974  independently  derived  equivalent  TSP  formulations  of  the  MTSP  [5], 
[58],  [65]  and  consequently  showed  that  the  M - salesmen  problem  is  no  more 
difficult  than  its  one-salesman  counterpart.  The  equivalence  is  obtained  by 
creating  m copies  of  the  origi.i,  each  connected  to  the  other  nodes  exactly 
as  was  the  original  origin  (with  the  same  distances).  The  ra  copies  are  not 
connected,  or  are  connected  by  arcs  each  with  distances  exceeding 

e”  I"  this  way,  an  optimal  single  salesman  tour  in  the  expanded 

i=l  j=l 

network  will  never  use  an  arc  connecting  two  copies  of  the  origin.  Then,  by 
coalescing  the  copies  back  into  a single  node,  the  single  salesman  tour  decom- 
poses into  M subtours  as  required  in  the  MTSP. 
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Vehicle  Dispatching 

The  vehicle  dispatching  problem  was  first  considered  by  Dantzig  and 
Ramser  [20]  who  developed  a heuristic  approach  using  linear  programming  ideas 
and  aggregation  of  nodes.  The  problem  is  to  obtain  a set  of  delivery  routes 
from  a central  depot  to  the  various  demand  points,  each  of  which  has  known 
requirements,  tc  minimize  the  total  distance  covered  by  the  entire  fleet. 

Vehicles  have  capacities  and  maximum  route  time  constraints.  All  vehicles 
start  and  finish  at  the  central  depot.  We  will  refer  to  the  following  formulation 
for  this  problem  as  the  generic  vehicle  routing  problem  (VRP) : 


Minimize 


subject  to 


" k " 


i=l  j=l 


k 
c . . 
ij 


n 

n 

NV 

1 

i=l 

Z 

j=l 

Z d, 

k=l  ’ 

k 

X . 

(3.1) 

n 

r 

i=l 

NV 

Z 

k=l 

k 

X . . 

ij 

= 1 

(j  - 2,  . 

n)  (3.2) 

n 

NV 

r 

Z 

Z 

IV 

X . . 

= 1 

(i  = 2,  . 

...  n)  (3.3) 

j=l 

k=l 

n 

Z 

i=l 

k 

X . 

n 

- Z 
j=l 

k 

h 

0 (k  = 1, 
P = 1, 

. . . , NV;  (3.4) 

. . . , n) 

n 

n 

Z 

i=l 

j 

2 X.  . ; 

=1 

< 

\ (k  = 1, 

...,  NV)(3.5) 

n 

Z 

i=l 

n 

Z 

j=l 

k k 

t . . X . 
XJ  i; 

■ T, 

] - ^ 

(k  = 1, 

...,  NV)(3.6) 

n 

Z 

j=2 

k 

"ij 

< 1 

(k  = 1, 

...,  NV)(3.7) 
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k 

I.  X.,  1 1 (k  = 1,  MV)  (3.8) 

i=2 

X e S (3.9) 


ij 


= 0 or  1 


for  all  i,  j,  k. 


(3.10) 


where 


n = number  of  nodes 
NV  = number  of  vehicles 
P|^  = capacity  of  vehicle  k 

= maximum  time  allowed  for  route  of  vehicle  k 
= demand  at  node  i (Q  = 0) 


t^  = time  required  for  vehicle  k to  deliver  or  collect  at  node  ? 

t . . = travel  time  for  vehicle  k from  node  i to  node  j (t^  = <») 
ij  ^ ii 


(t^=0) 


d^j  ® distance  from  node  i to  node  j 
k 


1 if  arc  (i,  j)  is  traversed  by  vehicle  k 
0 otherwise 

NV 


X = matrix  with  components  x..  = E x..  , specifying  connections 


ij 


regardless  of  vehicle  type. 


k=l 


Equation  (3.1)  states  that  total  distance  is  to  be  minimized.  Alternatively, 
we  could  minimize  costs  by  replacing  d^j^^  by  the  cost  coefficient  c^^  which 
depends  upon  the  vehicle  type.  Equations  (3.2)  and  (3.3)  ensure  that  each 
demand  node  is  served  by  one  vehicle  and  only  one  vehicle.  Route  continuity 
is  represented  by  equations  (3.4),  i.2.,  if  a vehicle  enters  a demand  node, 
it  must  exit  from  that  node.  Equations  (3.5)  are  the  vehicle  capacity  con- 
straints; similarly,  equations  (3.6)  are  the  total  elapsed  route  time  con- 
straints. For  instance,  a newspaper  delivery  truck  may  be  restricted  -om 
spending  more  than  one  hour  on  a tour  in  order  that  the  maximum  time  interval 
from  press  to  street  be  made  as  short  as  possible.  Equations  (3.7)  and  (3.8) 
make  certain  that  vehicle  availability  is  not  exceeded.  Finally,  the  subtour- 


breaking  constraints  (3.9)  can  be  any  of  the  equations  specified  previously. 

Since  (3.2)  and  (3.4)  imply  (3.3),  and  (3.4)  and  (3.7)  imply  (3.8),  from  now 

on  we  consider  the  generic  model  to  include  (3.1)  - (3.10)  excluding  (3.3) 

and  (3.8),  which  are  redundant.  We  assume  that  max  Q.  < min  P,  . 

l<i<„  W<k<NV 

That  is,  the  demand  at  each  node  does  not  exceed  the  capacity  of  any  truck. 
Observe  that  when  vehicle  capacity  constraints  (3.5)  and  route  time  constraints 
(3.6)  are  nonbinding,  and  can  be  ignored,  this  model  reduces  to  a multiple 
traveling  salesmen  problem  by  eliminating  constraints  (3.4)  and  substituting 
NV  ^ 

X..  = S X..  in  the  objective  function  and  remaining  constraints. 
k=l 


In  our  generic  model  we  have  assumed  that  when  a demand  node  is  serviced, 

its  requirements  are  satisfied.  In  other  words,  one  visit  is  sufficient.  A 

mixed  integer  programming  heterogeneous  fleet  problem  formulation  was  given 

in  1967  by  Garvin  [25]  in  which  this  assumption  is  relaxed.  The  number  of 

variables  in  his  formulation  is  much  greater  than  in  the  previous  model  where 

the  structure  is  related  more  closely  to  the  formulation  of  the  fundamental 

Traveling  Salesman  Problem.  Balinski  and  Quandt  [4]  provide  another  formulation 

in  terms  of  a set  cove  ing  problem. 

" k 

Observe  that  since  Z x . is  1 or  0 depending  upon  whether  or  not  the 

j=2 

th  k 

k vehicle  is  used,  a fixed  acquisition  cost  f.  4 . can  be  added  to  the 

j=l  ^ 


model  when  it  is  formulated  with  an  objective  function  to  investigate  tradeoffs 
between  routing  and  acquisition  costs. 

We  can,  without  difficulty,  generalize  the  generic  model  to  the  multi- 


commodity case  where  several  different  types  of  products  must  be  routed 
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simultaneously  over  a network  in  order  to  satisfy  demands  at  delivery  points 
for  the  various  products  [29]. 

Furthermore,  we  can  incorporate  timing  restrictions  into  the  vehicle 
dispatching  model.  If  we  define  a^  as  the  arrival  time  at  node  j then 
restrictions  on  delivery  deadlines  and  earliest  delivery  times  a^  can  be 
represented  by  the  following  nonlinear  equations: 


- - / . k . k 

aj  - E a 


(j  = 1 n) 


aj^  = 0 


a.  < a.  < a, 
-J  - J - j 


(j  = 2,  ...,  n). 


Alternatively,  the  above  nonlinear  constraints  can  be  replaced  by  the 
linear  constraints 


a . > (a.  + t^  + t^. ) 

J - i i 


(1  - Xjj)  T 


(1  - X..)T 


for  all  i,  j , k, 


where  T = max  T,  . Wlien  x , =0,  these  constraints  are  redundant.  When 

1 < k<  NV  ^ 


X,,  = 1,  they  determine  a.  in  terms  of  the  arrival  time  a.  at  the  node  i pre- 
i3  J 1 

ceeding  node  j on  a tour,  the  delivery  time  t^  at  node  i,  and  the  travel  time 
t^^  between  nodes  i and  j.  Of  course,  these  constraints  add  considerably 
to  the  size  of  the  model  (3.1)  - (3.10). 
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Multi-Depot  Vehicle  Routing 

The  integer  programming  formulation  of  the  vehicle  routing  problem  is 
altered  in  minor  ways  to  incorporate  multiple  depots.  Letting  nodes  1,  2, 

M denote  the  depots,  we  obtain  the  formulation  by  changing  the  index  in 
constraints  (3.2)  and  (3.3)  to  (j  = M+1,  ...,  n) , by  changing  constraints 
(3.7)  and  (3.8)  to 

M n . 

Z Z • 1 1 (k  = 1,  ....  MV)  (3.7') 

i=l  j=-MH  ^ 

M n 

Z Z X 1 1 (k  = 1 NV)  (3.8'), 

p=l  i=M+l 

and  by  redefining  our  previous  choices  for  the  subtour  breaking  set  S to  be: 

(!')  S = {(x..):  Z Z X..  ^ 1 ■ for  every  proper  subset 

ieQ  j^Q  Q of  V containing  nodes 

1,  2,  ...,  m}j 

(2')  S = {(x..):  Z Z X..  < Iq|-1  for  every  nonempty  sub- 

ieQ  jcQ  set  Q of  {^H-l,  Mf2,  ..., 

(3')  S = {(x^. ^):  y^.  - y^  + n x_. . < n-1  for  M-fl  £ i j £ n for 

some  real  numbers  y . } . 

^ 1 


Concluding  Comments  on  Formulations 

Hopefully,  the  formulations  discussed  so  far  provide  some  kind  of 
unified  basis  for  viewing  vehicle  routing  problems.  The  formulations  indicate, 
for  example,  that  Lagrangian  approaches  similar  to  those  used  by  Held  and 
Karp  [32],  [33]  can  be  applied  to  these  problems  by  dualizing  with  respect  to 
constraints  (3.2)  - (3.8).  Procedures  of  this  nature  might  prove  useful  when 


I 
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combined  with  good  heuristics.  In  particular,  since  the  heuristic  procedures 
provide  feasible  solutions  to  these  problems  and  the  Lagrangian  approaches 
provide  lower  bounds  to  the  minimum  distance  solution,  the  heuristic  and 
dual  approaches  together  bracket  the  optimal  value  for  the  objective  function. 
This  bracketing  might  conceivably  help  to  evaluate  the  effectiveness  of  the 
heuristic  approaches,  or  might  be  useful  in  obtaining  optimal  solutions  to  the 
problems  via  branch  and  bound  methods. 

III.  HEURISTIC  SOLUTION  TECHNIQUES  FOR  SINGLE  DEPOT  PROBLEMS 

Proposed  techniques  for  solving  vehicle  routing  problems  have  fallen 
into  two  classes  - those  which  solve  the  problem  optimally  by  branch  and 
bound  techniques,  and  those  which  solve  the  problem  heuristically . Since 
the  optimal  algorithms  have  been  viable  only  for  very  small  problems,  we 
concentrate  on  heuristic  algorithms.  Christofides  claims  that  the  largest 
vehicle  routing  problem  of  any  complexity  that  has  been  solved  exactly 
involved  only  23  customers  [17].  The  three  vehicle  routing  heuristic  methods 
which  we  will  discuss  (Clarke  and  Wright  [18],  Tyagi  [72],  and  Gillett  and 
Miller  [28])  have  been  used  for  problems  with  up  to  1000  customers. 

The  Clarke-Wright  algorithm  is  an  "exchange"  algorithm  in  the  sense 
that  at  each  step  one  set  of  tours  is  exchanged  for  a better  set  of  tours. 
Initially,  we  suppose  that  every  two  demand  points  i and  j are  supplied 
individually  from  two  vehicles  (refer  to  Figure  I.  below). 


Figure  I.  Initial  Setup. 


W 
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Now  if  instead  of  two  vehicles,  we  used  only  one,  then 
a savings  in  travel  distance  of 

= d,.  + d,.  - d,.  (see  Figure  IT.  below). 

} 1 Ij  ij  “ 


we  would  experience 

+ d, . + d. .) 

Ij  iJ 


Figure  II.  Nodes  i and  j have  been  linked. 


For  every  possible  pair  of  demand  points  i and  j there  is  a corresponding 
savings  . We  order  these  savings  from  greatest  to  least  and  starting  from 
the  top  of  the  list  we  link  nodes  i and  j where  s^^  represents  the  current 
maximum  savings  unless  the  problem  constraints  are  violated.  Christofides 
and  Eilon  found  from  10  small  test  problems  that  tours  produced  from  the 
"savings"  method  averaged  only  3.2  percent  longer  than  the  optimal  tours  [15]. 

Tyagi  [72]  presents  a method  which  groups  demand  points  in  the  following 
very  straightforward  fashion.  Starting  with  node  2 (nooa  1 is  the  central 
depot)  we  find  its  nearest  neighbor,  say  node  k,  subject  to  the  restriction 
that  Q.,  + Q.  £ C (C  is  capacity  of  the  vehicle  being  routed).  VJc  next 
find  the  nearest  neighbor  to  node  k,  say  node  j,  such  that  + Q.  + Q.  C 
and  continue  until  adding  a nearest  neighbor  wi71  result  in  a tour  exceeding 
either  vehicle  capacity  or  maximum  tour  length.  Rules  of  thumb  are  specified 
to  minimize  the  frequency  with  which  a group  will  consist  of  only  one  delivery 
point,  especially,  in  the  case  where  the  Uelivery  is  small  or  the  distance 
from  the  central  depot  to  this  point  is  more  than  half  the  distance  from  the 
farthest  point  to  the  central  depot.  Having  grouped  the  delivery  points  into 
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m tours,  the  vehicle  dispatching  problem  red.ices  to  m Traveling  Salesman 
Problems,  one  for  each  tour. 

In  a recent  paper,  Gillett  and  Miller  [28]  introduce  an  efficient 
buildup  algorithm  for  handling  up  to  about  250  nodes.  Rectangular  coordinates 
for  each  demand  point  are  required,  from  which  we  may  calculate  polar  coordinates. 
We  select  a "seed"  node  randomly.  With  the  central  depot  as  the  pivot,  we 
start  sweeping  (clockwise  or  counterclockwise)  the  ray  from  the  central 
depot  to  the  seed.  Demand  nodes  are  added  to  a route  as  they  are  swept.  If 
the  polar  coordinate  indicating  angle  is  ordered  for  the  demand  points  from 
smallest  to  largest  (with  seed's  angle  0)  we  enlarge  routes  as  we  Increase 
the  angle  until  capacity  restricts  us  from  enlarging  a route  by  including  an 
additional  demand  node.  This  denuind  point  becomes  the  seed  for  the  following 
route.  Once  we  have  the  routes  we  can  apply  TSP  algorithms  such  as  the 
Lin-Kernighan  heuristic  to  improve  tours  and  obtain  significantly  better 
results.  In  addition,  we  can  vary  the  seed  and  select  the  best  solution. 

The  Clarke  and  Wright  algorithm  overcomes  a major  deficiency  of  the  other 
two  algorithms  in  that  demand  points  farther  av/ay  from  the  central  depot  are 
considered  early  for  linking.  This  property  means  that  few  of  these  "problem" 
nodes  are  left  to  be  grouped  at  the  end  when  there  are  virtually  no  degrees 
of  freedom.  On  the  other  hand,  Gillett  and  Miller  form  non-overlapping  tours, 
unlike  the  other  two  algorithms.  The  Tyagi  algorithm  is  fast  and  the  easiest 
to  program. 

The  literature  contains  computational  experience  relating  to  the  Clarke 
and  Wright  [18],  and  the  Gillett  and  Miller  [28]  algorithms.  On  a 50-customer 
problem  (which  will  be  mentioned  later)  the  former  algorithm  had  computation 
time  of  6 seconds,  the  latter,  120  seconds.  However,  the  second  objective 
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value  was  slightly  lower.  A 250-location  problem  with  an  average  of  10 
locations  per  route  was  solved  in  just  under  10  minutes  of  IBM  360/67  time 
using  the  Gillett  and  Miller  algorithm.  The  Tyagi  algorithm  has  been 
programmed  by  Klincewicz  [Al].  The  same  50-node  problem  was  solved  in  1 second 
of  execution  time  yielding  a rather  high  objective  value.  Several  experiments 
with  the  Tyagi  algorithm  have  confirmed  the  poor  objective  performance  of 
this  approach.  Yellow's  modification  of  the  Clarke-Wr ight  procedure  appears 
computationally  to  be  the  most  powerful  vehicle  routing  method  (based  on 
computational  experience  mentioned  in  the  literature);  problems  of  200  nodes 
have  been  solved  in  less  than  a minute  and  a problem  with  1000  nodes  was 
solved  in  five  minutes  on  an  IBM  360/50  [80].  Webb,  in  applying  a similar 
sequential  savings  approach,  reports  having  solved  AOO-customer  problems  in 
less  than  6 seconds  on  a GDC  6600  [75].  Several  questions  arise  since 
implementation  is  not  discussed  in  the  paper  at  all: 

(i)  Are  input  and  output  operations  included  in  this  figure? 

(ii)  Can  the  program  handle  different  types  of  vehicles? 

(iii)  Is  there  a maximum  number  of  drop  points  per  route  in  the  program, 
or  a maximum  travel  time  for  the  vehicles? 

(iv)  Do  the  computations  involve  integer  or  real  arithmetic? 

(v)  Although  the  GDC  machine  is  faster  than  the  IBM  machine,  how 
does  one  explain  the  extraordinary  difference  in  reported 
running  times  between  Yellow's  and  Webb's  papers? 

(vi)  How  is  the  data  stored? 
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IV.  A MODIFIED  CLARKE-WRIGHT  ALGORITHM  WHICH  IS  VERY  FAST 


Probably  the  most  popular  of  the  heuristic  solution  techniques  discussed 
in  the  last  section  is  the  Clarke-Wright  "savings"  method  which  IBM  has 
programmed  as  VSPX  [36],  a flexible  computer  code  to  handle  complex  routing 
problems. 

At  each  step  in  the  Clarke-Wright  algorithm  we  seek  the  greatest  positive 
savings  s^^  subject  to  the  following  restrictions: 

(i)  nodes  i and  j are  not  already  on  the  same  route; 

(ii)  neither  i nor  j are  Interior  to  an  existing  tour; 

(iii)  vehicle  availability  is  not  exceeded; 

(iv)  vehicle  capacity  is  not  exceeded; 

(v)  maximum  number  of  drops  (or  maximum  route  time)  is  not  exceeded. 

In  many  applications,  where  the  delivery  or  pickip  time  is  a sizeable  portion 
of  total  route  time,  it  makes  sense  to  limit  the  number  of  drops  rather  than  the 
route  time  since  travel  times  are  so  difficult  to  estimate.  Nodes  i and  j, 
then,  are  linked  together  to  form  a new  route,  and  the  procedure  is  repeated 
until  no  further  savings  are  possible. 

There  have  been  modifications  made  to  the  basic  approach,  most  notably 
those  suggested  by  Beltrami  and  Bodin  [7]  and  by  Yellow  [80].  In  this  paper, 
we  emphasize  data  structures  and  list  processing  and  discuss  a new  implementation 
of  the  Clarke-Wright  algorithm  which  is  motivated  by  (i)  optimality  considerations, 
(ii)  storage  considerations,  and  (iii)  sorting  considerations  and  program 
running  time.  We  modify  the  basic  Clarke-Wright  algorithm  in  three  ways: 

(1)  by  using  a route  shape  parameter  Y to  define  a modified  savings 

s . . = d,  . + d,  . - Y d , . 

ij  lx  Ij  ij 


-17- 


and  fioding  the  best  route  structure  obtained  as  the  parameter 
is  varied; 

(2)  hy  considerine  savings  only  between  nodes  that  are  "close" 
to  each  other; 

(3)  by  storing  savings  s^.  in  a heap  structure  to  reduce  comparison 
operations  and  ease  aicess. 

Route  Shape  Parameter 

The  route  shape  parameter  was  introduced  by  Yellow  [80],  and  is  an  out- 
growth of  an  algorithm  developed  by  Gaskell  [26].  Wlien  the  parameter  y Increases 
from  zero,  greater  emphasis  is  placed  on  the  distance  between  points  i and  j 
rather  than  their  position  relative  to  the  central  depot.  The  search  for  the 
best  route  structure  as  Y is  varied  provides  heuristic  solutions  which  are 
closer  to  the  optimal  solution  than  otherwise  obtained  via  the  traditional 
algorithm  where  Y = !•  ^or  example,  in  applying  the  modified  algorithm  to 
the  50  node  problem  at  Christofides  and  Eilon  [15],  a solution  of  total 
distance  577  was  obtained  with  a route  shape  parameter  of  1.3,  whereas  the 
traditional  Clarke-Wright  solution  is  585.  Perhaps  more  importantly,  the 
search  over  y promotes  flexibility.  In  most  applications  there  are  "unstated" 
goals  and/or  constraints.  By  varying  the  route  shape  parameter  we  can  produce 
several  "sufficiently  good"  tentative  solutions  from  which  a final  solution  may 
be  chosen  on  the  basis  of  these  unstated  considerations. 

Storage 

The  Clarke-Wright  algorithm  was  designed  initially  to  handle  a matrix 
of  real  inputs,  distances  and  savings.  One  might  argue  that  these  variables 
can  be  rounded  to  integers  to  simplify  computations.  However,  we  felt  that 
precision  was  important.  In  dealing  with  a 600  node  problem  a matrix  form 
would  require  360,000  storage  locations  for  these  inputs.  This  assumes  an 
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undirected  network  with  symmetric  distances  in  order  that  the  d^^  terms 
can  fill  above  the  diagonal  in  the  matrix,  and  the  s^^  terms  can  fill  below. 

Of  course,  if  we  use  only  the  traditional  Clarke-Wright  algorithm  we  need  not 
store  the  distances  at  all.  At  least  some  distances  must  be  kept  in  memory, 
however,  in  order  to  vary  the  route  shape  parameter  (in  particular,  the  distances 
from  the  origin  to  all  other  nodes  must  be  stored). 

Rather  than  consider  all  pairwise  linkings,  as  in  a matrix  approach, 
we  can  become  more  selective  and  work  with  the  linkings  of  greatest  interest 
and  convenience.  We  have  been  investigating  a Newspaper  distribution  problem 
of  nearly  600  nodes;  our  computer  system  cannot  set  aside  360,000  internal 
storage  locations  (more  than  1.5  million  bytes).  And  if  it  were  possible,  it 
would  be  inefficient  to  calculate  savings  which  could  never  be  realized.  Most 
of  the  potential  linkings  are  infeasible  from  a practical  standpoint  and 
our  modified  Clarke-Wright  program  reflects  this  observation.  The  topology 
of  the  network  is  stored  in  ladder  representation  form.  In  general,  for  each 

arc  we  record  its  origin  node  and  its  destination  node,  and  its  length. 

2 

The  approach  requires  3 A rather  than  n storage  locations  where  n is  the  total 
number  of  nodes  and  A the  number  of  arcs  under  consideration.  For  undirected 
networks  we  can  order  the  arcs  lexicographically  and  cut  ladder  representation 
storage  in  half. 

Given  a grid  of  width  WIDTH  and  height  HEIGHT,  and  the  x and  y coordinates 

2 

of  the  nodes  of  the  transportation  network,  the  grid  is  divided  into  DIV 
rectangles  in  such  a way  that  each  demand  node  is  contained  in  a rectangle  of 
width  WIDTH/DIV  and  height  HEIGHT/DIV.  Node  I has  box  coordinates  BX(I)  and 


BY  (I).  The  set  of  arcs  we  will  oe  working  with  includes  all  arcs  between  the 
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central  depot  (node  1)  and  other  nodes  (demand  nodes),  and  all  arcs  between 

nc les  which  are  no  further  apart  than  one  box.  In  other  \;ords,  if 

|BX(I)  - BX(J)1  > 1 or  |BY(I)  - BY(J)1  > 1 the  arc  (I,  J)  is  ignored. 

The  value  of  DIV  influencet  the  accuracy  of  the  heuristic  algorithm  and  should 
be  altered  according  to  the  problem.  If  the  number  of  demand  points,  n-1,  is 
small  DIV  should  be  small;  if  n is  large  DIV  should  be  large.  The  smaller 
the  parameter  DIV,  for  a particular  problem,  the  larger  is  the  number  of 
arcs  to  be  considered. 

Heap  Structure 

At  each  step  of  the  algorithm,  we  must  determine  the  maximum  savings.  This 

comparison  of  savings  can  be  accomplished  quickly  and  conveniently  by  partially 

ordering  the  data  in  a heap  structure  and  updating  the  structure  at  each 

step  after  altering  the  routes.  More  precisely,  the  savings  s^^  as  denoted 

by  Sj^,  ...,  s^  are  arranged  in  a binary  tree  with  k levels  called  a heap. 

The  essential  property  of  a heap  is  that  s.  > s_,  and  s.  > s_.,,.  Below  is 

1—  2i  1—  2x+l 

a heap  for  k = 3 (s^^  = 25,  s^  = 21,  s^  = 16,  and  so  on). 


Figure  III. 


If  the  list  does  not  completely  fill  the  last  level  of  the  tree,  then 
we  can  add  positions  in  the  last  level  with  entries  of  - °°.  Clearly, 
corresponds  to  the  maximum  savings  possible  of  those  savings  rnder  consideration. 
First,  the  subroutine  STHEAP  arranges  Sj^,  ...,  into  a heap  using 


r 
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che  ideas  in  Williams  [77]  and  Floyd  [23],  and  produces  a vector  JX«  points 
to  where  positions  from  the  pre-heap  listing  of  savings  are  in  the  current 
heap.  For  example,  if  initially  an  element  was  in  position  64  and  currently 
it  occupies  position  32  then  JX(64)  = 32.  STHEAP  builds  a heap  in  0 (log2m) 
comparisons  and  interchanges  in  the  worst  case.  Now,  suppose  that  s^^  corres- 
ponds to  arc  (i,  j)  and  that  neither  i nor  j become  interior  points  on  a route 
upon  linking  nodes  i and  j.  Setting  s^^  to  zero  we,  in  effect,  remove  s, 
from  the  heap  since  only  positive  savings  are  considered.  Actually,  we 
bubble  Sj^  down  the  heap  until  it  finds  its  new  home.  A new  heap  can  be 
constructed  in  this  case  with  remarkable  ease  by  subroutine  OTHEAP  since  s^^ 
is  always  made  smaller  and  must  therefore  move  down  the  binary  tree.  If  Sj^ 
has  already  been  moved  to  position  i then  the  new  entry  s^  is  compared  with 
its  sons  only  (82^  and  82^^^^  until  s^  ^ max  (s2^,  ®2i-H^‘ 

If,  on  the  other  hand,  node  i (or  j)  becomes  an  interior  point,  then  we 
no  longer  try  linking  node  i (or  3)  with  any  other  nodes  in  the  network. 
Conceptually,  we  can  cross  all  entries  s^  off  the  heap,  where  s^  is  the  savings 
incurred  by  linking  an  interior  node  to  another  node.  We  acconplish  this 
task  by  (i)  setting  each  such  entry  to  zero  and  (ii)  having  subroutine 
OTHEAP  reconstruct  a heap  structure.  OTHEAP  is  called  as  many  times  as  there 
are  adjacent  nodes  to  the  interior  node.  The  matrix  ADJA  (I,  J)  indicates 

the  position  on  the  pre-heap  listing  where  the  potential  savings  incurred 

th 

from  linking  node  I with  its  J adjacent  neighbor  can  be  obtained.  The 
vector  J>‘  enables  us  to  determine  the  current  position  of  this  entry  on  the 
heap;  this  vector  is  being  updated  throughout  the  computer  progrr*" 

The  first  entry  of  the  heap  always  represents  the  most  promising  link 
to  add  to  the  current  routes,  and  by  eliminating  infeasible  linkings  quickly 


:) 
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by  OTHEAP,  tlip  algorithm  proceeds  very  rapidly  until  no  adaicional  savings 
can  be  realized.  OTHEAP  exploits  the  simple  observation  that  when  an  entry 
is  changed  it  is  always  decreased  to  zero.  The  stei->s  of  the  OTHEAP  algorithm 
are  given  below. 

Basic  OTHEAP  Procedure 

Step  0.  NN  is  the  length  of  the  heap.  IC  is  the  positio  . of  the  entry 
which  has  been  decreased  to  zero. 

Step  1.  II  TC 

COPY  S(II). 

Step  2.  J II  + II. 

Step  3.  If  node  II  has  no  son  on  the  heap  then  go  to  Step  7. 

Step  4.  If  max  {S(J),  S(J+1)}  = S(J+1)  then  J *■  J+1. 

Step  5.  If  S(J)  £ COPY  go  to  Step  7. 

Step  6.  S(II)  ^ S(J) 

II  J 

Go  to  Step  2. 

Step  7.  S(II)  ^ COPY 

Stop. 

Computational  Results 

The  computer  program  has  been  programmed  in  the  WATFIv  version  of  bUKTKAN. 
Computation.  1 studies  have  beer  especially  encouraging.  A 50  node  problem 
(mentioned  earlier)  was  solved  on  the  IBM  370/168  at  M.I.T.  using  less  than 
one  second  total  execution  time,  including  all  input  and  output  operations. 
Nodes  were  read  in  by  coordinates,  .so  all  distances  were  calculated  within  the 
program.  Realistic  data  for  a newspaper  distribution  problem  was  gathered 
for  us  by  a local  newspaper.  The  evening  edition  has  a city  circulation  of 
about  100,000  with  nearly  600  demand  points  where  bundles  of  newspapers  are 
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delivered.  The  total  execution  time  was  20  seconds  and  the  routes  produced 
were  considered  from  reasonable  to  very  interesting  by  the  circulation 
department  involved.  We  believe  our  program  can  be  of  considerable  value 
as  a tool  in  helping  to  alleviate  newspaper  distribution  problems.  It  is 
of  interest  that  a major  portion  of  the  20  seconds  involved  the  determination 
of  node  adjacencies  and  distances. 

Below  we  report  sensitivity  analyses  for  the  newspaper  problem  with 
respect  to  the  following  parameters: 

(i)  box  sizes  (Table  I); 

(ii)  vehicle  capacities  (Table  II); 

(iii)  maximum  number  of  drops  on  a tour  (Table  III) ; 

(iv)  route  shape  parameter  (Table  IV) . 

In  each  case,  only  one  parameter  is  varied  at  a time.  The  tendencies  exhibited 
in  Tables  I,  II,  and  III  agree  with  our  intuition,  whereas  Table  IV  is  a bit 
more  interesting.  The  best  choice  for  route  shape  parameter  appears  to  be 
truly  problem  dependent  (in  this  example,  y ==  1*0  and  Y = are  strong 
candidates,  taking  both  objective  value  and  number  of  tours  into  account). 
Several  parameters  might  be  varied  simultaneously  in  future  parametric  studies. 
This  parametric  analysis  affords  us  the  opportunity  to  choose  the  route 
structure  most  appealing  in  light  of  various  (often  conflicting)  objectives. 
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Value 
of  DIV 

Total  distance 
traveled 

Number  of 
routes 

30 

306.19 

22 

25 

272.55 

17 

24 

259.67 

15 

23 

267.11 

17 

22 

261.83 

16 

Table  I.  Testing  box  sizes. 


Vehicle 

capacity 


Total  distance 
traveled 


Number  of 
routes 


• 

3000 

266.20 

18 

3500 

262.32 

17 

AOOO 

262.32 

17 

4500 

262.32 

17 

' 

5000 

262.32 

17 

Table  II.  Testing  vehicle  capacities  (the  capacity  of  the  smaller  of  two 
types  of  vehicles  was  varied) . 


Maximum 

number  of  drops 

Total  distance 
traveled 

Number  of 
routes 

60 

262.32 

17 

70 

262.65 

16 

80 

260.36 

15 

90 

259.86 

15 

100 

259.67 



15 

Table  III.  Testing  the  maximum  number  of  drops  on  a tour. 


Value 
of  Y 

Total  distance 
traveled 

Number  of 
routes 

.4 

263.09 

15 

.6 

264.72 

16 

.8 

266.03 

17 

1.0 

262.32 

17 

1.2 

277.46 

18 

1.4 

.279.77 

19 

1.6 

276.99 

20 

Table  IV.  Testing  the  route  shape  parameter. 


1 

. V.  HEURISTIC  SOLUTION  TECHNIQUES  FOR  MULTI-DEPOT  PROBLEMS 

f 

I - V'rtiereas  the  VRP  has  been  studied  widely,  the  multi-depot  problem  has 

I attracted  less  attention.  The  relevant  literature  is  represented  by  only  a 

few  papers.  The  exact  methods  for  the  single  depot  VRP  can  be  extended  to  the 
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general  ease,  but  as  in  the  case  for  the  single  depot  VRP,  because  of  storage  and 
computation  time,  only  small  problems  can  be  dealt  with  currently.  Problems 
arising  in  practice  are  beyond  the  scope  of  these  optimal  algorithms.  Three 
relatively  successful  heuristic  approaches  have  been  developed  for  the  multi- 
depot problem. 

A first  class  of  heuristics  generates  one  solution  arbitrarily  and  proceeds 
to  improve  the  solution  by  exchanging  nodes  one  at  a time  between  routes  until 
no  further  improvement  is  possible.  Some  typical  examples  are  Wren  and  Holliday 
[79]  and  Cassidy  and  Bennett  [13]. 

The  Wren  and  Holliday  program  has  been  run  on  Gaskell's  ten  sample  cases 
and  the  authors  report  results  comparable  with  other  methods  with  respect 
to  computation  time  and  accuracy.  No  tests  involving  more  than  100  nodes  are 
reported.  Cassidy  and  Bennett  report  a successful  case  study  using  a similar 
approach. 

The  two  other  approaches  are  extensions  of  heuristics  developed  for  the 
single  depot  VRP.  The  Gillett  and  Johnson  algorithm  [27]  is  an  extension 
of  the  Gillett  and  Miller  "sweep"  algorithm  discussed  previously.  The  method 
solves  the  multi-depot  problem  in  two  stages.  First,  locations  are  assigned 
to  depots.  Then,  several  single  depot  VRP's  are  solved.  Each  stage  is  treated 
separately. 

Initially,  all  locations  are  in  an  unassigned  state.  For  any  node  i,  let 
t'  (i)  be  the  closest  depot  to  i,  and 
t''(i)  be  the  second  closest  depot  to  i. 

For  each  node  i,  the  ratio 


r(i) 


1,  t (i)  1,  t (i) 


t 
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is  computed  and  all  nodes  are  ranked  by  Increasing  values  of  ^(i).  The 
arrangement  determines  the  ordei  in  which  the  nodes  are  assigned  to  a depot: 
those  which  are  relatively  close  to  a depot  are  assigned  first  and  those 
lying  midway  between  two  depots  are  considered  last.  After  a certain  number 
of  nodes  have  been  assigned  from  the  list  of  r(i),  a small  cluster  is 
formed  around  every  depot.  Locations  i such  that  the  ratio  r(i)  is  close 
to  1,  are  assigned  more  carefully. 

If  two  nodes  j and  k are  already  assigned  to  a depot  t,  inserting  i 
between  j and  k on  a route  linked  to  t creates  an  additional  distance  equal  to 


d.  . + d., 
Ji  ik 


which  represents  a part  of  the  total  distance  (or  cost).  In  other  words,  the 
algorithm  assigns  location  i to  a depot  t by  inserting  i between  tv;o  nodes 
already  assigned  to  depot  t,  in  the  least  costly  manner. 

The  Sweep  Algorithm  is  used  to  construct  and  sequence  routes  in  the 
cluster  about  each  depot  independently.  A number  of  refinements  are  made  to 
improve  the  current  solution. 

Gillett  and  Johnson  [27]  provide  a list  of  problems  and  the  corresponding 

solutions  given  by  their  algorithm.  We  will  make  use  of  this  list  in  order 

to  evaluate  our  multi-depot  algorithm.  Let 

d..  = distance  between  demand  nodes  i and  i,  and 

d^  = distance  between  node  i and  depot  k. 


The  algorithm  starts  with  the  initial  solution  consisting  of  servicing 


each  node  exclusively  by  one  route  from  the  closest  depot.  The  total  distance 


N 

of  all  routes  is  then  D = Z 2min 

i=l  k 


The  method  successively 


links  pairs 


(* 
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of  cities  in  order  to  decrease  the  total  distance  traveled.  One  basic  rule 
is  assumed  in  the  algorithm:  the  first  assignment  of  nodes  to  the  nearest 

terminal  is  temporary  but  once  two  or  more  nodes  have  been  assigned  to  a 
common  route  from  a termina? , the  nodes  are  not  reassigned  to  another 
terminal. 

At  each  step,  the  choice  of  linking  a pair  of  nodes  i and  j on  a route 
from  terminal  k is  made  in  terms  of  two  criteria: 

(i)  a savings  when  linking  i and  j at  k, 

(ii)  a penalty  for  not  doing  so. 

Nodes  i and  j can  be  linked  only  if  no  constraints  are  violated. 

The  computation  of  savings  is  not  as  straightforward  as  in  the  case  of 
a single  depot.  The  savings  s^^  are  associated  to  every  combination  of 
demand  nodes  i and  j and  depot  k,  and  represent  the  decrease  in  total  distance 
traveled  obtained  when  linking  i and  j at  k.  The  formula  is  given  by 


= d"?  + - d.. 

1 J ij 


(4) 


where 


d'^ 

1 


r 


2 min 


d/ 

V.  ^ 


if  i has  not  yet  been  given 
a permanent  assignment 

otherwise. 


An  illustration  of  formula  (4)  is  given  in  Nguyen  [54]. 

The  savings  s^^  are  computed  for  i,  j = 1,  ...,  N (i  ^ j)  and 

k = 1,  ...,  M at  each  step.  They  can  be  stored  in  M matrices,  each  N by  N. 

Tillman  and  Cain  [68]  introduce  a penalty  correction  factor  to  the 
Clarke-Wright  procedure  as  follows.  Suppose  that  nodes  i and  j are  linked 
on  a route  from  depot  k at  a particular  iteration  with  savings  s, . . 


If,  instead. 
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J 


I 


I 


nodes  i and  u were  linked  on  a route  from  depot  m at  a later  iteration,  the 

m Ic 

savings  would  be  s^^.  Assuming  that  s^^  is  selected  to  give  maximum  savings 


jj  ^ 

out  of  node  i,  s.  < s.,  and  0. 

’ iu  — ij  xu 


km 

s^^  - s^^  measures  the  opportunity 


cost  of  using  this  later  link  addition.  Similarly,  if  at  a later  iteration, 
node  i were  linked  to  node  v on  a route  from  depot  m,  then  0 ™ = s.^  - s°. 

would  be  the  corresponding  opportunity  cost.  Using  these  observations, 

I. 

Tillman  and  Cain  define  a penalty  p^.  for  not  selecting  link  (1,  j)  on  a 

••■J 

route  from  depot  k at  the  current  iteration  by: 

p^j  = min  {0^™|  all  (u,m),  u=l,  ...»  N;  m=l,  ...,  M,  except  (j,k)) 

+ min  (0^™|  all  (v , m),  v=l,  ...,  N;  m=l,  ...,  M,  except  (i,  k)}. 

In  order  to  incorporate  the  concepts  of  savings  and  penalties  in  a simple 
way,  Tillman  and  Cain  suggest  the  expression 


c k ^ . Q k 

'ij  ■ “=ij  " ^^ij 


where  a and  3 are  two  selected  (or  varied)  positive  weights.  At  each  step, 
the  link  (i,  j)  at  k is  chosen  which  maximizes  f^^  and  which  yields  a feasible 
solution  (with  regard  to  capacity  and  maximum  route  time  restrictions) . 

The  idea  of  considering  a penalty  function  can  obviously  be  applied  to 
the  single  depot  VRP  as  an  improvement  on  Clarke  and  Wright's  original  algorithm. 
The  latter  has  the  drawback  of  proceeding  very  myopically  and  the  consideration 
of  a penalty  function  might  help  much  in  overcoming  this  drawback.  The  penalty 
calculations,  however,  impose  additional  computational  burdens. 
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VI.  ALGORITHM  I FOR  THE  MULTI-DEPOT  VRP 

The  algorithm  proposed  here  is  based  on  the  savings  method,  and  is 
motivated  by  a desire  to  find  an  algorithm  capable  of  handling  large  size 
problems.  Algoritlira  I uses  Tillman  and  Cain's  approach  for  computing  savings 
but  excludes  the  idea  of  a penalty  function.  Our  main  contribution,  at  this 
stage,  is  a method  of  manipulating  data  which  allows  fast  computation  and 
minimizes  storage  requirements. 

Algorithm  I extends  Clarke  and  Wright's  algorithm  to  the  case  of  many 
depots.  The  idea  is  to  start  with  an  initial  solution  consisting  of  each 
node  being  served  exclusively  by  a route  from  the  closest  depot.  The  algorithm 
proceeds  by  selecting  at  each  step  the  largest  savings  s^^^ , computed  as  in 
Tillman  and  Cain  [68],  and  i and  j are  linked  on  a route  served  from  depot  k. 
Once  a link  has  been  created  it  is  never  removed  in  later  steps,  and  the 
assignment  of  nodes  i and  j to  depot  k is  also  final. 

In  most  codes  written  for  the  VRP  using  the  savings  approach,  the  problem 
of  storing  the  list  of  savings  otren  has  set  the  limit  on  the  size  of  problems 
tractable.  In  the  multi-depot  VRP,  to  each  depot  corresponds  a matrix  of 
savings  and  the  total  number  of  storage  locations  needed  grows  quickly 
(even  though  the  number  of  depots  is  generally  small  for  transportation 
networks).  In  order  to  cope  with  the  problem  of  storage,  we  have  used  the 
following  approach. 

!iC 

First,  we  note  that  the  savings  associated  with  the  linkage  of 

nodes  i and  j at  k is  symmetric  in  i and  j . We  can  exploit  this  symmetry  by 
storing  savings  in  rectangular  rather  than  square  matrix  form.  We  perform 
the  operation  illustrated  in  Figures  IV  and  V.  Ni  is  defined  to  be  N/2  if 
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N is  even  and  (N+l)/2  if  N is  odd  (N  being  the  number  of  demand  nodes  in 
the  problem) . We  replace  elements  of  triangle  I by  elements  of  triangle 
II,  row  by  row,  starting  from  the  last  row  of  triangle  II,  and  in  reverse 
order.  We  obtain  a rectangular  matrix  Sj^  of  dimension  (N-1)  x N1  which 
contains  all  the  necessary  data.  Examples  for  N even  and  odd  are  given  in 
Figure  VI. 


redundant 

information 


Figure  IV.  Storing  savings  in  a square  matrix. 


Each  element  (I,  J) 
which  is  stored  in  location 

For  N even,  Cl,  J)  = 

For  N odd,  Sj^  (I,  J)  = 


of  the  matrix  Sj^  corresponds  to  a savings  s^^ 
(I,  J)  of  Sj^  according  to  the  following  rules: 


N-I+1,  N-J+1 


if  I ^ J 
if  I < J; 


\ I+l,  J if  I > J 

/ 

/ N-I+1,  N-J+1  if  I < J < (N-D/2 

\ 0 if  I < J = (N+D/2. 


For  large  problems,  a further  economy  is  achieved  by  transforming  the 
into  half-word  integers,  caution  being  taken  against  overflow.  Round-off 
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triangle  1 


triangle  II 


Figure  V.  Storing  savings  in  rectangular  form. 


N. 

"7,5 

0 
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"5,3 

®6,1 

"6,2 

®6,3 

"6,4 

®7,i 

"7,2 

"7,3 

"7,4 

s 

2,1 

\8,7 

®8,6 

"8,5 

®3,1 

"3,2 

\"7,6 

"7,5 

G 

"4,1 

"4,2 

"4,3\ 

"6,5 

• 

• 

• 

s\ 

5,4  ^ 

• 

• 

• 

"6,4 

• 

• 

• 

"7,4 

®8,1 

®8,2 

"8,3 

"8,4 

Figure  VI.  Examples  of  savings  matrices. 
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errors  can  be  minimized  by  first  multiplying  all  real  distances  by  a factor 
of  100. 

In  the  computer  program,  we  define  ROWMAX(I)  as  the  largest  element  in 
all  the  rows  I of  all  the  M matrices  Sj^,  k=l,  . . . , M.  At  each  step  of  the 
algorithm,  the  maximal  savings  is  obtained  by  searching  for  the  maximum 
among  the  ROWIAX(I) , 1=1,  ...,  N-1.  This  can  be  readily  accomplished  via 
heap  structures,  as  we  have  noted  previously.  The  main  advantage  of  using 
the  auxiliary  variable  ROWMAX  and  searching  for  the  maximal  savings  among 
ROTOIAX(I),  1=1,  ....  N-1  instead  of  searching  directly  through  the  whole 
list  of  savings  is  an  economy  of  computation  time.  At  each  iteration  of  the 
algorithm,  we  have  to  update  some  rows  and  columns  of  the  savings  matrix. 
Typically,  this  operation  involves  changes  only  for  a few  elements  of  the 
vector  RO^'JMAX.  Updating  the  vector  ROWMAX  and  then  selecting  its  largest 
element  involves  much  less  comparisons  than  searching  through  the  whole 
list  of  new  savings. 

Each  node  is  represented  by  its  rectangular  coordinates.  Four  pointers 
are  associated  with  each  node  I:  L(I),  RT(I),  TERM(I),  and  SQ(I).  We 

define  these  below: 

( 2 if  I is  the  only  node  on  a route 

L(I)  = j 1 if  I is  an  endpoint  on  a route 

I 0 if  I is  an  interior  point  on  a route 

RT(I)  is  the  number  of  the  route  which  contains  I 

TERM(I)  is  the  terminal  to  which  I is  assigned 
SQ(I)  is  the  load  of  the  truck  which  serves  I. 

The  algorithm  starts  with  a solution  consisting  of  each  node  being  exclusively 
served  by  a route  from  the  closest  depot  and  at  each  subsequent  step  we  link  node 


I 


i and  j at  terminal  k in  order  to  realize  the  largest  possible  savings. 


Steps  of  Algorithm  I 


Step  0. 


Step  1. 


L(I)  = 2 
RT(I)  = I 
SQ(I)  = Q(l) 


for  I ® 1,  . . . , N. 


For  i=l,  N,  find  d^*^.  Set  TERM(I)  c and  MIND(I)  = d^*^ 

where  MIND(I)  is  the  distance  from  terminal  t to  demand  node  1. 


Step  2. 


Step  3. 


Compute  savings: 

f 0 if  SQ(I)  + SQ(J)  > CAP 

/ 2 MIND(I)  + 2 MIND(J)  - d/  - d.*^  - d if  SQ(I)  + SQ(J)  < CAP. 

L ^ J ij 

Determine  ROWIAX(I)  for  I^lf  N-1,  and  for  each  I,  determine 

the  indices  KMAX(I)  and  JMAX(l)  such  that  ROWMAX(I)  = 

W(i)  <'■ 

Compute  OPT  = Max  {ROWMAX(I) | 1=1,  N-1}.  lOPT  is  the  index 

of  OPT  (i.e.,  OPT  = ROmAX(IOPT)) . If  OjPT  = 0,  go  to  Step  . 

If  OPT  > 0,  trace  back  to  the  savings  s,,  which  corresponds  to 
OPT  as  follows:  ^ 

Let  JOPT  = JMAX(IOPT)  and  KOPT  = KMAX(IOPT). 

If  lOPT  > JOPT  then  -T  i = lOPT  + 1, 

/ j = JOPT,  and 

KOPT 

Here,  OPT  = 


’ IOPT+1,  JOPT' 

If  lOPT  < JOPT  then  \ i = N - lOPT  + 1, 

1 j = N - JOPT  + 1,  and 


„KOPT 


Here,  OPT  - _ jQpx  + i* 


Nodes  1 and  j are  linked  at  terminal  k. 


Step  5.  Update  pointers: 

TERM(I)  = TERM(J)  = k. 

L(I)  L(I)  - 1. 

L(J)  L(J)  - 1. 

For  all  nodes  w e RT(I)  or  RT(J),  RT(w)  = RT(I)  and 

SQ(w)  = SQ(I)  + SQ(J). 

DIST  <-  DIST  + d , . . 

ij 
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Perform  Step  6 for  h = i and  then  h = j . 

Step  6.  If  L(h)  = 1: 
in 

set  0 for  all  rngf^k  and  all  w; 

set  s|^^  sj^^  + 2 d|^  - 2 MIND(h)  for  all  w (wj^i) 

if  SQ(h)  + SQ(w)  £ CAP; 

set  S,  -t-  0 if  SQ(h)  + SQ(w)  > CAP. 
hw 


Step  7. 
Step  8. 


If  L(h)  = 0: 

set  S™  0 for  all  m and  all  w. 
hw 

Update  ROWIAX(I)  for  1=1,  N-1. 

Go  to  Step  A. 

Compute  total  distance 

N-1 

DIST  DIST  + Z Ui)  . d and  then  stop. 

i=l  ^ 


Modifications  and  Improvements 

The  version  of  Algorithm  I we  have  described  requires  the  computation  of 

all  savings  for  which  Q(i)  + Q(j)  £ CAP  at  the  initialization  step.  The 

N 

number  of  computations  here  can  be  close  to  M(^).  However,  it  is  not  necessary 
to  search  through  all  possible  linkages  of  i to  other  nodes.  It  is  sufficient 
to  consider  only  those  nodes  which  are  close  neighbors  of  i.  We  set  up 
a grid  here  as  described  in  a pievious  section. 

In  order  to  assess  the  ;,conomy  of  computation  brought  about  by  the  grid, 
we  have  tried  the  idea  on  a problem  with  50  nodes  and  A depots:  without 
using  a grid.  Algorithm  I took  2.0A  seconds  total  CPU  time  and  with  a grid, 
it  took  1.38  seconds  while  giving  the  same  solution. 

As  with  the  single  depot  VRP  a route  shape  parameter  line  search  has 
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been  studied.  In  general  there  seems  to  be  no  way  a priori  to  determine 
a best  value  of  y for  any  given  problem.  The  algorithm  is  fast  enough, 
however,  so  that  ve  can  try  different  values  of  Y and  then  selec*"  the  best 
solution  produced.  We  have  tried  the  idea  on  a 50-node,  4-depot  problem; 
results  are  given  in  Table  V. 


Value 

Solution  produced  by  algorithm 

of  a 

Total  distance 

Number  of 

traveled 

routes 

.2 

584.63 

6 

.4 

559.89 

6 

.6 

512.09 

6 

.8 

509.77 

6 

1.0 

509.30 

7 

1.2 

508.07 

7 

1.4 

518.10 

7 

1.6 

536.40 

11 

1.8 

560.40 

13 

2.0 

587.28 

15 

Table  V.  Testing  route  shape  parameter  for  multi-depot  problems. 

The  advantage  of  this  approach  is  that  we  are  given  a few  alternative.*? 
to  chose  from  and  depending  o'<  the  objective  function,  one  solution  or 
another  can  be  selected.  For  example,  in  the  above  illustration,  the  sole 
criterion  of  minimal  distance  traveled  would  determine  the  choice  of 
Y = 1.2;  but  the  consideration  of  distance  traveled  combined  with  the  number 
of  routes  might  set  the  value  Y ” -8  as  the  best  choice. 

Another  positive  feature  of  the  method  is  that  after  drawing  the 
different  routes  produced  for  different  values  of  Y,  it  is  usually  possible 
to  combine  some  of  them  manually,  simply  by  examining  the  various  alternatives. 
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( 


to  product  an  even  better  solution.  For  example,  we  worked  a few  minutes  on 
the  solution  obtained  with  Y = 1.0,  and  obtained  a solution  with  6 routes 
and  a total  distance  traveled  of  487  units.  We  will  describe  a computer 
code  which  performs  a similar  refining  operation  later  in  this  paper. 

The  above  idea  is  of  course  a simplistic  application  of  the  interactive 
approach  to  VRP  solving.  The  appeal  of  the  method  is  its  possibility  of 
combining  extensive  and  long  computations  done  with  a computer  with  the 
intuition  and  judgement  of  the  human  mind.  Krolak  et.  al.  [44]  have 
conducted  research  in  this  direction. 


Some  Computaticnal  Results 


We  have  run  a program  implementing  the  above  ideas  for  some  problems 
which  were  solved  by  Gillett  and  Johnson  [27].  A version  of  Algorithm  I 
which  is  especially  well-suited  for  problems  with  an  even  number  of  depots 
was  used.  Instead  of  storing  the  savings  in  matrices  ..f  dimension  (N-1)  x Nl, 
we  store  them  in  M/2  (M  is  even)  square  matrices  each  N x N;  the  upper  half 
of  the  k*”^’  matrix  corresponds  -o  savings  associated  with  depot  2k  - 1,  and 
the  lower  half  corresponds  to  savings  associated  with  depot  2k.  Clearly, 
an  odd  number  M of  depots  creates  no  additional  snags  unless  M is  large. 

This  version  is  slightly  faster  than  the  original  Algorithm  I since  less 
manipulation  of  indices  is  involved. 
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Figure  VII.  Savings  matrix  for  M even. 

The  code  we  have  written  using  the  above  storage  scheme  is  compiled 
in  0.18  seconds  on  an  IBM  370/168  and  takes  8400  bytes  of  core.  A problem 
of  200  nodes  and  4 depots  necessitates  200k  of  core  memory  (using  halfword 
integers  for  savings),  300  nodes  and  4 depots  take  400k,  200  nodes  and  10 
depots  take  420k.  We  consider  these  to  be  medium-sized  problems.  Results 
for  some  sample  problems  are  summarized  in  Table  VI. 


VII.  ALGORITHM  II:  EXTENSION  OF  ALGORITHM  I FOR  LARGE  PROBLEMS 

With  larger  problems,  for  example  with  1000  nodes,  an  approach  is  to 
divide  the  problem  into  as  many  subproblems  as  there  are  depots  and  to 
solve  each  subproblem  separately.  Basically,  there  are  two  steps  involved 
in  a muJti-depot  VRP:  first,  nodes  have  to  be  allocated  to  depots;  then 

routes  are  built  which  link  nodes  assigned  to  the  same  depot.  Ideally,  it 
is  more  efficient  to  deal  with  the  two  steps  simultaneously,  as  Algorithm  I 
does,  but  the  method  is  no  longer  computationally  tractable  when  the  number 
of  nodes  becomes  too  large. 
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The  idea  behind  Algoritlim  II  is  that  if  a given  node  i is  much  closer 

to  one  depot  than  any  other  depot,  i will  be  served  from  its  closest  depot. 

U'lien  node  i is  equidistant  from  several  depots,  the  assignment  of  i to  a 

depot  becomes  more  difficult.  For  every  node  i,  we  determine  the  closest 

k k 

depot  k,  and  the  second  closest  depot  k_.  If  the  ratio  r = d.  ^/d.  2 is 
i L i i i 

less  than  a certain  chosen  parameter  6 (0  ^ 6 £ 1) , wa  assign  i to  k^^; 
in  the  case  r^  ^ 6 we  say  that  node  i is  a border  node.  The  use  of  the  ratio 
r^  appears  in  Gillett  and  Johnson  [27]  also. 

Clearly,  if  6 = 0,  all  nodes  are  declared  border  nodes  and  if  6 = 1, 
all  nodes  are  assigned  to  their  closest  depot.  For  a given  problem,  we 
can  fix  the  number  of  border  nodes  as  we  wish,  by  varying  6. 

The  method  proceeds  as  follows.  In  the  first  step,  we  ignore  the  non- 
border nodes  and  Algorithm!  I is  applied  to  the  set  of  border  nodes.  The 
algorithm  allocates  the  border  nodes  to  depots  and  simultaneously  builds 
segments  of  routes  connecting  these  nodes.  >\t  the  end  of  this  first  step, 
all  nodes  of  our  problem  are  assigned  to  some  depot  and  all  border  nodes 
are  linked  on  some  routes.  The  solution  to  the  VRP  is  produced  depot  by 
depoc  using  single  depot  VRP  techniques.  The  segments  of  routes  which  are 
built  on  border  nodes  are  extended  to  the  remaining  nodes. 

It  is  felt  that  the  efficiency  of  the  method  depends  on  how  many  border 
nodes  Algorithm  I can  handle.  One  approach  which  allows  a maximal  number 
of  border  nodes  to  be  considered  involves  the  ordering  of  nodes  by  decreasing 
ratios  r^.  Alternatively,  one  can  experiment  with  a list  of  increasing 
values  of  6.  As  6 increases  from  0 to  1 the  number  of  border  nodes 
decreases.  The  method  described  here  has  been  implemented  with  real  data 
taken  from  the  distribution  information  of  a local  nev;spaper.  The  problem, 
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with  almost  600  nodes  and  2 depots,  ran  in  under  55  seconds. 


VIII.  POST-PROCESSOR 


In  a previous  section  we  mentioned  the  possibility  of  improving  the 
solution  obtained  by  Algorithm  I by  modifying  the  routes  that  it  produced. 
In  this  section,  we  discuss  a computerized  procedure  to  perform  this  task. 

We  suppose  that  initially  we  are  given  a solution  to  the  VRP.  An 
arbitrary  orientation  is  assigned  to  each  route  so  that  for  each  node  i, 
we  can  define  pr(i)  to  be  the  node  or  depot  which  precedes  i on  its  route 
and  fl(i)  to  be  the  node  or  depot  which  follows  i on  its  route. 

If  a node  j is  inserted  between  nodes  i and  fl(i),  the  reduction  in 
distance  traveled  can  be  computed  as: 


*^pr(j),  j ^ fl(j)  ^ '^i,  fl(i)  "^prO),  fl(j)  ^i,  j ^ j , fl(i) 


To  any  pair  of  nodes  i and  j , we  can  consider  the  savings  u . . corresponding 

ij 

to  the  insertion  of  node  j between  i and  fl(i).  In  addition,  if  node  i 

is  the  first  node  served  by  a route,  then  the  savings  v^^  associated  with 

the  insertion  of  node  j between  i and  pr(i)  has  to  be  taken  into  account. 

If  the  number  of  routes  is  r,  then  the  possible  savings  can  be  stored  in 

a rectangular  matrix  W of  dimension  (N  + r)  x N.  The  first  N rows  correspond 

to  savings  u..  and  the  last  r rows  refer  to  the  savings  v . . . We  note  that 

in  general  w. . / w. . . 

® ij  ji 

The  post-processor  searches  for  the  largest  element  of  W which  is 


feasible,  and  performs  the  indicated  insertion.  The  operation  is  then  repeated 
until  no  further  reduction  in  total  distance  traveled  is  possible.  The  details 
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are  very  much  like  those  of  Algorithm  I and  are  not  repeated  here. 

The  post-processor  uses  a savings  approach  which  differs  from  the 
Clarke  and  Wright  algorithm  in  several  respects  [54].  In  general,  updating 
the  savings  matrix  is  more  cumbersome  as  many  more  savings  might  change  in 
value.  Typically,  in  the  Clarke  and  Wright  algorithm,  iterations  get  shorter 
and  shorter  as  the  procedure  progresses  because  the  number  of  nodes  to 
consider  decreases,  whereas  in  the  post-processor,  they  remain  about  the 
same  in  length. 

With  a code  for  implementing  the  ideas  discussed,  we  have  tried  the 
approach  with  problems  reported  in  Table  VI.  Results  are  summarized  in 
Table  VII. 

Our  method  has  produced  solutions  comparable  with  the  ones  reported 
by  Gillett  and  Johnson  while  using  much  less  computation  time.  In  an 
exceptional  case  the  Gillett  and  Johnson  algorithm  required  227.67  seconds 
of  CPU  time;  our  solution  was  slightly  higher  while  using  only  2-3%  of  the 
time. 

IX.  CONCLUSION 

This  paper  has  presented  the  vehicle  routing  problem  and  the  multi-depot 
VRP,  discussed  various  integer  programming  formulations,  studied  several 
heuristic  approaches,  and  introduced  very  efficient  algorithms  which  have 
been  implemented  successfully.  The  indication  is  that  the  suggested  procedures 
are  efficient  and  can  be  used  as  effective  decision-making  tools  for  large 


scale  vehicle  routing  problems  encountered  in  practice. 
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